home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / randist / tdist.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-05-14  |  2.1 KB  |  81 lines
  1. /* randist/tdist.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <math.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_sf_gamma.h>
  24. #include <gsl/gsl_rng.h>
  25. #include <gsl/gsl_randist.h>
  26.  
  27. /* The t-distribution has the form
  28.  
  29.    p(x) dx = (Gamma((nu + 1)/2)/(sqrt(pi nu) Gamma(nu/2))
  30.    * (1 + (x^2)/nu)^-((nu + 1)/2) dx
  31.  
  32.    The method used here is the one described in Knuth */
  33.  
  34. double
  35. gsl_ran_tdist (const gsl_rng * r, const double nu)
  36. {
  37.   if (nu <= 2)
  38.     {
  39.       double Y1 = gsl_ran_ugaussian (r);
  40.       double Y2 = gsl_ran_chisq (r, nu);
  41.  
  42.       double t = Y1 / sqrt (Y2 / nu);
  43.  
  44.       return t;
  45.     }
  46.   else
  47.     {
  48.       double Y1, Y2, Z, t;
  49.       do
  50.     {
  51.       Y1 = gsl_ran_ugaussian (r);
  52.       Y2 = gsl_ran_exponential (r, 1 / (nu/2 - 1));
  53.  
  54.       Z = Y1 * Y1 / (nu - 2);
  55.     }
  56.       while (1 - Z < 0 || exp (-Y2 - Z) > (1 - Z));
  57.  
  58.       /* Note that there is a typo in Knuth's formula, the line below
  59.      is taken from the original paper of Marsaglia, Mathematics of
  60.      Computation, 34 (1980), p 234-256 */
  61.  
  62.       t = Y1 / sqrt ((1 - 2 / nu) * (1 - Z));
  63.       return t;
  64.     }
  65. }
  66.  
  67. double
  68. gsl_ran_tdist_pdf (const double x, const double nu)
  69. {
  70.   double p;
  71.  
  72.   double lg1 = gsl_sf_lngamma (nu / 2);
  73.   double lg2 = gsl_sf_lngamma ((nu + 1) / 2);
  74.  
  75.   p = ((exp (lg2 - lg1) / sqrt (M_PI * nu)) 
  76.        * pow ((1 + x * x / nu), -(nu + 1) / 2));
  77.   return p;
  78. }
  79.  
  80.  
  81.